In this script I’ll analyze the genetic clustering done with plink.
Preliminaries
Calling libraries
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(GGally)##
## Attaching package: 'GGally'## The following object is masked from 'package:dplyr':
##
## nasalibrary(ggpubr)## Loading required package: magrittrReading data
setwd('..')
path <- getwd()
setwd(paste(path, "/Results/GenClustering", sep = ""))
pca.eigenval <- read.csv("Clus.eigenval", sep = "", header = F)
pca.eigenvec <- read.csv("Clus.eigenvec", sep = "", header = F)Data cleaning
Naming data columns
PCnames <- sprintf("PC%s",seq(1:50))
colnames(pca.eigenvec) <- c("FID", "IID", PCnames)Creating new column to plot hapmap pops against our samples.
pca.eigenvec$database <- pca.eigenvec$FID
pca.eigenvec$database[!(pca.eigenvec$FID == "ASW" | pca.eigenvec$FID == "CEU" | pca.eigenvec$FID == "CHB" | pca.eigenvec$FID == "CHD" |
pca.eigenvec$FID == "GIH" | pca.eigenvec$FID == "JPT" | pca.eigenvec$FID == "LWK" | pca.eigenvec$FID == "MEX" |
pca.eigenvec$FID == "MKK" | pca.eigenvec$FID == "TSI" | pca.eigenvec$FID == "YRI")] <- "ADAPT"Because plink output is a normalized PC, clustering did not achieved very good results. Now we de-normalized the PCA by multiplying de eigenvectors by the standard deviation.
pca.eigenvec[,3:52] <- t(c(as.matrix(pca.eigenval)) * t(as.matrix(pca.eigenvec[,3:52])))I’ll remove samples that contain the word genome (surely those are not in the face files), the word ver (those are controls from the are from the axiom arrays), the word gDNA, as well as those with the name A_103. I’ve seen that some of those are outliers in some PCs
removed <- filter(pca.eigenvec, IID == "A_103" | grepl("^genome", IID) |
grepl("^ver", IID) | grepl("^gDNA", IID)) %>%
select(contains("ID"))
removed## FID IID
## 1 1001 genome_001
## 2 1002 genome_002
## 3 1003 genome_003
## 4 1004 genome_005
## 5 1005 genome_006
## 6 1006 genome_007
## 7 1007 genome_009
## 8 1008 genome_010
## 9 1009 genome_013
## 10 1010 genome_016
## 11 1010 genome_11002_1
## 12 1011 genome_018
## 13 1011 genome_Brian_Schimmoller_Full_20131206130740
## 14 1012 genome_019
## 15 1012 genome_Jennifer_K_Wagner_Full_20130924125453
## 16 1013 genome_020
## 17 1013 genome_JKW
## 18 1014 genome_022
## 19 1014 genome_William_Mattern_Full_20140619082357
## 20 1015 genome_024
## 21 1016 genome_026
## 22 1017 genome_027
## 23 1018 genome_028
## 24 1019 genome_033
## 25 1020 genome_038
## 26 1021 genome_039
## 27 1022 genome_042
## 28 1023 genome_047
## 29 1024 genome_049
## 30 1025 genome_051
## 31 1026 genome_053
## 32 1027 genome_056
## 33 1028 genome_061
## 34 1029 genome_062
## 35 1030 genome_064
## 36 1031 genome_074
## 37 1032 genome_076
## 38 1033 genome_079
## 39 1034 genome_081
## 40 1035 genome_083
## 41 1036 genome_094
## 42 1037 genome_096
## 43 1038 genome_097
## 44 1039 genome_101
## 45 1040 genome_104
## 46 1041 genome_105
## 47 1042 genome_107
## 48 1043 genome_109
## 49 1044 genome_112
## 50 1045 genome_114
## 51 1046 genome_117
## 52 1047 genome_118
## 53 1048 genome_119
## 54 1049 genome_120
## 55 1050 genome_121
## 56 1051 genome_122
## 57 1052 genome_123
## 58 1053 genome_127
## 59 1054 genome_132
## 60 1055 genome_137
## 61 1056 genome_139
## 62 1057 genome_141
## 63 1058 genome_142
## 64 1059 genome_144
## 65 1060 genome_147
## 66 1061 genome_154
## 67 1062 genome_158
## 68 1063 genome_173
## 69 1064 genome_186
## 70 1065 genome_195
## 71 1066 genome_196
## 72 1067 genome_197
## 73 1068 genome_200
## 74 1069 genome_201
## 75 1070 genome_202
## 76 1071 genome_203
## 77 1072 genome_204
## 78 1073 genome_205
## 79 1074 genome_206
## 80 1075 genome_207
## 81 1076 genome_208
## 82 1077 genome_209
## 83 1078 genome_210
## 84 1079 genome_211
## 85 1080 genome_212
## 86 1081 genome_213
## 87 1082 genome_214
## 88 1083 genome_215
## 89 1084 genome_216
## 90 1085 genome_217
## 91 1086 genome_218
## 92 1087 genome_219
## 93 1088 genome_220
## 94 1089 genome_221
## 95 1090 genome_222
## 96 1091 genome_223
## 97 1092 genome_224
## 98 1093 genome_225
## 99 1094 genome_226
## 100 1095 genome_227
## 101 1096 A_103
## 102 1096 genome_229
## 103 1097 genome_230
## 104 1098 genome_232
## 105 1099 genome_233
## 106 10100 genome_234
## 107 10101 genome_235
## 108 10102 genome_236
## 109 10103 genome_237
## 110 10104 genome_238
## 111 10105 genome_239
## 112 10106 genome_240
## 113 10107 genome_241
## 114 10108 genome_242
## 115 10109 genome_243
## 116 10110 genome_244
## 117 10111 genome_245
## 118 10112 genome_246
## 119 10113 genome_247
## 120 10114 genome_248
## 121 10115 genome_249
## 122 10116 genome_250
## 123 10117 genome_938
## 124 10118 A_103
## 125 10119 A_103
## 126 10120 A_103
## 127 10121 A_103
## 128 10122 A_103
## 129 10123 A_103
## 130 10124 A_103
## 131 10125 A_103
## 132 10126 A_103
## 133 AffyAxiom gDNA_103
## 134 AffyAxiom ver10
## 135 AffyAxiom ver11
## 136 AffyAxiom ver12
## 137 AffyAxiom ver13
## 138 AffyAxiom ver14
## 139 AffyAxiom ver15
## 140 AffyAxiom ver16
## 141 AffyAxiom ver17
## 142 AffyAxiom ver18
## 143 AffyAxiom ver19
## 144 AffyAxiom ver20
## 145 AffyAxiom ver21
## 146 AffyAxiom ver22
## 147 AffyAxiom ver23
## 148 AffyAxiom ver24
## 149 AffyAxiom ver25
## 150 AffyAxiom ver26
## 151 AffyAxiom ver27
## 152 AffyAxiom ver28
## 153 AffyAxiom ver29
## 154 AffyAxiom ver30
## 155 AffyAxiom ver31
## 156 AffyAxiom ver32
## 157 AffyAxiom ver33
## 158 AffyAxiom ver34
## 159 AffyAxiom ver35
## 160 AffyAxiom ver37
## 161 AffyAxiom ver39
## 162 AffyAxiom ver41
## 163 AffyAxiom ver42
## 164 AffyAxiom ver43
## 165 AffyAxiom ver44
## 166 AffyAxiom ver45
## 167 AffyAxiom ver46
## 168 AffyAxiom ver47
## 169 AffyAxiom ver48
## 170 AffyAxiom ver49
## 171 AffyAxiom ver50
## 172 AffyAxiom ver51
## 173 AffyAxiom ver52
## 174 AffyAxiom ver53
## 175 AffyAxiom ver54
## 176 AffyAxiom ver55
## 177 AffyAxiom ver56
## 178 AffyAxiom ver57
## 179 AffyAxiom ver58
## 180 AffyAxiom ver59
## 181 AffyAxiom ver60
## 182 AffyAxiom ver61
## 183 AffyAxiom ver62
## 184 AffyAxiom ver63
## 185 AffyAxiom ver64
## 186 AffyAxiom ver65
## 187 AffyAxiom ver66
## 188 AffyAxiom ver67
## 189 AffyAxiom ver68
## 190 AffyAxiom ver69
## 191 AffyAxiom ver70
## 192 AffyAxiom ver71
## 193 AffyAxiom ver72
## 194 AffyAxiom ver73
## 195 AffyAxiom ver74
## 196 AffyAxiom ver75
## 197 AffyAxiom ver76
## 198 AffyAxiom ver77
## 199 AffyAxiom ver78
## 200 AffyAxiom ver79
## 201 AffyAxiom ver80
## 202 AffyAxiom ver81
## 203 AffyAxiom ver82
## 204 AffyAxiom ver83
## 205 AffyAxiom ver84
## 206 rando genome_Amie_Atiyeh_20090419204251
## 207 rando genome_Brett_Patterson_Full_20100313124936
## 208 rando genome_Charles_Voss_Full_20100313013353
## 209 rando genome_Daris_Wells_Full_20100317110224
## 210 rando genome_Emily_Voss_Full_20100313073416
## 211 rando genome_Henry_Louis_Gates_20090402070340
## 212 rando genome_Jeannette_Henney_Full_20100228133658pca.eigenvec <- filter(pca.eigenvec, IID != "A_103" & !grepl("^genome", IID) &
!grepl("^ver", IID) & !grepl("^gDNA", IID))I removed 212 samples.
We also have some outliers in PC5, and PC10. Probably there are also in some other lower PCs
p1 <- ggplot(pca.eigenvec, aes(PC5)) + geom_histogram()
p2 <- ggplot(pca.eigenvec, aes(PC10)) + geom_histogram()
ggarrange(p1, p2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
pca.eigenvec %>% filter(PC5 > 1) %>% select(contains("ID"))## [1] FID IID
## <0 rows> (or 0-length row.names)pca.eigenvec %>% filter(PC10 < -0.5) %>% select(contains("ID"))## FID IID
## 1 MKK NA21399
## 2 MKK NA21402
## 3 MKK NA21405pca.eigenvec <- pca.eigenvec %>% filter(PC5 < 1 & PC10 > -0.5)Descriptive and exploratory analyses
Some summary stats
summary(pca.eigenvec[,c(1,2,53)])## FID IID database
## ADAPT :2536 61556 : 2 ADAPT :5133
## AffyAxiom : 920 61572 : 2 MKK : 140
## cv : 673 61634 : 2 YRI : 113
## SanAntonio: 225 62022 : 2 CEU : 112
## MKK : 140 62531 : 2 LWK : 90
## v3 : 116 62577 : 2 GIH : 88
## (Other) :1508 (Other):6106 (Other): 442PCA screeplot
barplot(c(as.matrix(pca.eigenval)), names.arg=PCnames,
main = "Variances",
xlab = "Principal Components",
ylab = "Eigenvalue",
col = "steelblue")
Looking at the HapMap populations:
| POP | Description |
|---|---|
| ASW | African ancestry in Southwest USA |
| CEU | Utah residents with Northern and Western European ancestry from the CEPH collection |
| CHB | Han Chinese in Beijing, China |
| CHD | Chinese in Metropolitan Denver, Colorado |
| GIH | Gujarati Indians in Houston, Texas |
| JPT | Japanese in Tokyo, Japan |
| LWK | Luhya in Webuye, Kenya |
| MXL | Mexican ancestry in Los Angeles, California |
| MKK | Maasai in Kinyawa, Kenya |
| TSI | Toscani in Italia |
| YRI | Yoruba in Ibadan, Nigeria |
plothap <- function(dat, x1, y1)
{
p1 <- ggplot(dat, aes_string(x = x1, y = y1, color = "database")) +
geom_point(alpha = 0.3, size = 2) + theme_pubr()
p1 <- ggpar(p1, palette = "jama")
return(p1)
}
p1 <- plothap(filter(pca.eigenvec, database != "ADAPT"), "PC1", "PC2")
p2 <- plothap(filter(pca.eigenvec, database != "ADAPT"), "PC3", "PC4")
p3 <- plothap(filter(pca.eigenvec, database != "ADAPT"), "PC5", "PC6")
p4 <- plothap(filter(pca.eigenvec, database != "ADAPT"), "PC7", "PC8")
p5 <- plothap(filter(pca.eigenvec, database != "ADAPT"), "PC9", "PC10")
p6 <- plothap(filter(pca.eigenvec, database != "ADAPT"), "PC11", "PC12")
ggarrange(p1, p2, p3, p4, p5, p6, nrow = 3, ncol = 2,
common.legend = TRUE, legend = "right",
align = "v")
meanshap <- filter(pca.eigenvec, database != "ADAPT") %>% group_by(database) %>%
select(contains("PC")) %>% summarise_all(funs(mean))## Adding missing grouping variables: `database`plothapmeans <- function(dat, x1, y1){
p1 <- ggplot(dat, aes_string(x = x1, y = y1, label = "database")) +
geom_label(aes(fill = database), colour = "white", fontface = "bold", alpha = 0.6) +
theme_pubr(legend = "none")
p1 <- ggpar(p1, palette = "jama")
return(p1)
}
p1 <- plothapmeans(meanshap, "PC1", "PC2")
p2 <- plothapmeans(meanshap, "PC3", "PC4")
p3 <- plothapmeans(meanshap, "PC5", "PC6")
p4 <- plothapmeans(meanshap, "PC7", "PC8")
p5 <- plothapmeans(meanshap, "PC9", "PC10")
p6 <- plothapmeans(meanshap, "PC11", "PC12")
ggarrange(p1, p2, p3, p4, p5, p6,
nrow = 3, ncol = 2,
align = "v")
If we plot our samples together with the HapMap, we see the following
plothapsamp <- function(dat1, dat2, x1, y1)
{
p1 <- ggplot(dat1, aes_string(x = x1, y = y1, label = "database")) +
geom_point(data = dat2,
aes_string(x = x1, y = y1), alpha = 0.3) +
geom_label(aes(fill = database), colour = "white", fontface = "bold", alpha = 0.4) +
theme_pubr(legend = "none")
p1 <- ggpar(p1, palette = "jama")
return(p1)
}
p1 <- plothapsamp(meanshap, filter(pca.eigenvec, database == "ADAPT"), "PC1", "PC2")
p2 <- plothapsamp(meanshap, filter(pca.eigenvec, database == "ADAPT"), "PC3", "PC4")
p3 <- plothapsamp(meanshap, filter(pca.eigenvec, database == "ADAPT"), "PC5", "PC6")
p4 <- plothapsamp(meanshap, filter(pca.eigenvec, database == "ADAPT"), "PC7", "PC8")
p5 <- plothapsamp(meanshap, filter(pca.eigenvec, database == "ADAPT"), "PC9", "PC10")
p6 <- plothapsamp(meanshap, filter(pca.eigenvec, database == "ADAPT"), "PC11", "PC12")
ggarrange(p1, p2, p3, p4, p5, p6,
nrow = 3, ncol = 2,
align = "v")
Looking at our samples with population designation
plothapgroup <- function(dat1, dat2, x1, y1)
{
p1 <- ggplot(dat1, aes_string(x = x1, y = y1, label = "database")) +
geom_point(data = dat2, aes_string(x = x1, y = y1, color = "FID"), alpha = 0.5) +
geom_label(aes(fill = database), colour = "white", fontface = "bold", alpha = 0.5) +
theme_pubr()
p1 <- ggpar(p1, palette = "jama")
return(p1)
}
eurhap <- filter(meanshap, database == "CEU" | database == "TSI")
eursamp <- filter(pca.eigenvec, FID == "Italian" | FID == "Polish" | FID == "Irish" |
FID == "Portuguese")
p1 <- plothapgroup(eurhap, eursamp, "PC1", "PC2")
p2 <- plothapgroup(eurhap, eursamp, "PC3", "PC4")
p3 <- plothapgroup(eurhap, eursamp, "PC5", "PC6")
p4 <- plothapgroup(eurhap, eursamp, "PC7", "PC8")
p5 <- plothapgroup(eurhap, eursamp, "PC9", "PC10")
p6 <- plothapgroup(eurhap, eursamp, "PC11", "PC12")
ggarrange(p1, p2, p3, p4, p5, p6,
nrow = 3, ncol = 2,
common.legend = TRUE, legend = "right",
align = "v")
Clustering
For an initial clustering attempt I’ll group the HapMap populations into continents, as shown in the table:
| POP | Continent |
|---|---|
| ASW | Africa |
| CEU | Europe |
| CHB | Asia |
| CHD | Asia |
| GIH | American |
| JPT | Asia |
| LWK | Africa |
| MXL | American |
| MKK | Africa |
| TSI | Europe |
| YRI | Africa |
pca.eigenvec$continent[(pca.eigenvec$database == "CEU" | pca.eigenvec$database == "TSI")] <- "Europe"
pca.eigenvec$continent[(pca.eigenvec$database == "ASW" | pca.eigenvec$database == "LWK" |
pca.eigenvec$database == "MKK" | pca.eigenvec$database == "YRI")] <- "Africa"
pca.eigenvec$continent[(pca.eigenvec$database == "CHB" | pca.eigenvec$database == "CHD" |
pca.eigenvec$database == "JPT")] <- "Asia"
pca.eigenvec$continent[(pca.eigenvec$database == "GIH" | pca.eigenvec$database == "MEX")] <- "America"If we look at scatter plots with the continent variable
plotcont <- function(dat, x1, y1)
{
p1 <- ggplot(dat, aes_string(x = x1, y = y1, color = "continent")) +
geom_point(alpha = 0.3, size = 2) + theme_pubr()
p1 <- ggpar(p1, palette = "jama")
return(p1)
}
p1 <- plotcont(filter(pca.eigenvec, database != "ADAPT"), "PC1", "PC2")
p2 <- plotcont(filter(pca.eigenvec, database != "ADAPT"), "PC3", "PC4")
p3 <- plotcont(filter(pca.eigenvec, database != "ADAPT"), "PC5", "PC6")
p4 <- plotcont(filter(pca.eigenvec, database != "ADAPT"), "PC7", "PC8")
ggarrange(p1, p2, p3, p4, nrow = 2, ncol = 2,
common.legend = TRUE, legend = "right",
align = "v")
I’ll use the HapMap samples as training, with a 5-fold cross-validation, and predict the values of our samples.
library(caret)## Loading required package: latticeset.seed(10)
train <- pca.eigenvec %>% filter(database != "ADAPT")
test <- pca.eigenvec %>% filter(database == "ADAPT")
# 5-fold repeated cross-validation
fitControl <- trainControl(method = "repeatedcv", number = 5, repeats = 5)
svm.fit <- train(continent ~ PC1 + PC2 + PC3 + PC4 + PC5, method = 'svmRadial',
trControl = fitControl, data = train)## Loading required package: kernlab##
## Attaching package: 'kernlab'## The following object is masked from 'package:ggplot2':
##
## alphasvm.fit## Support Vector Machines with Radial Basis Function Kernel
##
## 985 samples
## 5 predictor
## 4 classes: 'Africa', 'America', 'Asia', 'Europe'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times)
## Summary of sample sizes: 788, 788, 787, 789, 788, 788, ...
## Resampling results across tuning parameters:
##
## C Accuracy Kappa
## 0.25 0.9979685 0.9971561
## 0.50 0.9979685 0.9971561
## 1.00 0.9979685 0.9971561
##
## Tuning parameter 'sigma' was held constant at a value of 0.8521666
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.8521666 and C = 0.25.test$continent <- predict(svm.fit, test)
summary(test$continent)## Africa America Asia Europe
## 1185 647 215 3086Scatter plots of the predicted continent groups
p1 <- plotcont(test, "PC1", "PC2")
p2 <- plotcont(test, "PC3", "PC4")
p3 <- plotcont(test, "PC5", "PC6")
p4 <- plotcont(test, "PC7", "PC8")
ggarrange(p1, p2, p3, p4, nrow = 2, ncol = 2,
common.legend = TRUE, legend = "right",
align = "v")
The next step will be to create a clustering of suthern versus northern european populations, using CEU and TSI as training sets. We’ll double check our clustering with the samples that we have population information (Irish, Polish, Italian, and Portuguese).
pca.eigenvec$Eur[pca.eigenvec$database == "CEU"] <- "NEruope"
pca.eigenvec$Eur[pca.eigenvec$database == "TSI"] <- "SEruope"
set.seed(10)
train.eur <- pca.eigenvec %>% filter(database == "CEU" | database == "TSI")
test.eur <- filter(test, continent == "Europe")
# 5-fold repeated cross-validation
fitControl <- trainControl(method = "repeatedcv", number = 5, repeats = 5)
svm.fit.eur <- train(Eur ~ PC1 + PC2 + PC3 + PC4 + PC5, method = 'svmRadial',
trControl = fitControl, data = train.eur)
svm.fit.eur## Support Vector Machines with Radial Basis Function Kernel
##
## 200 samples
## 5 predictor
## 2 classes: 'NEruope', 'SEruope'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times)
## Summary of sample sizes: 160, 160, 160, 160, 160, 161, ...
## Resampling results across tuning parameters:
##
## C Accuracy Kappa
## 0.25 0.9160019 0.8304345
## 0.50 0.9170044 0.8323780
## 1.00 0.9129787 0.8240120
##
## Tuning parameter 'sigma' was held constant at a value of 0.2075051
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.2075051 and C = 0.5.test.eur$continent <- predict(svm.fit.eur, test.eur)
summary(test.eur$continent)## NEruope SEruope
## 1902 1184p1 <- plotcont(test.eur, "PC1", "PC2")
p2 <- plotcont(test.eur, "PC3", "PC4")
p3 <- plotcont(test.eur, "PC5", "PC6")
p4 <- plotcont(test.eur, "PC7", "PC8")
ggarrange(p1, p2, p3, p4, nrow = 2, ncol = 2,
common.legend = TRUE, legend = "right",
align = "v")
levels(test$continent) <- c("Africa", "America", "Asia", "Europe", "NEruope", "SEruope")
test[test$continent == 'Europe', 54] <- test.eur$continent
test$Eur <- NULL
test <- droplevels(test)Combining datasets
setwd('..')
path <- getwd()
setwd(paste(path, "/Results/PCA", sep = ""))
getwd()## [1] "C:/Users/tomas/Box Sync/Research/FacialSD/Results/PCA"#Read data
coeffs <- read.csv("coeffs.csv")colnames(coeffs)[1] <- "IID"
total <- merge(test, coeffs, by = "IID")
FacePCA <- total[,59:145]
GenPCA <- total[,3:52]
Covariates <- total[,c(1,2,55:58,54)]
PCnames <- sprintf("PC%s",seq(1:87))
colnames(FacePCA) <- PCnames
PCnames <- sprintf("PC%s",seq(1:50))
colnames(GenPCA) <- PCnames
MergedDat <- list(Covariates, FacePCA, GenPCA)
setwd('..')
path <- getwd()
setwd(paste(path, "/Results/MergedData", sep = ""))
save(MergedDat, file = "MergedDat.RData")